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ABSTRACT 

A two-dimensional, compressible Navier-Stokes 
equations with a k - e turbulence model are solved 
numerically to simulate the flows of compressible free 
shear layers. The appropriate form of k and i equations 
for compressible flows are discussed. Sarkar s modelling 
is adopted to simulate the compressibility effects in the 
k and e equations. The numerical results show that the 
the spreading rate of the shear layers decreases with 
increasing convective Mach number. In addition, fa- 
vorable comparison was found between the calculated 
results and Goebel and Dutton’s experimental data. 


INTRODUCTION 

Recent national interest in trans-atmospheric ve- 
hicle has rekindled the hypersonic research. For this 
vehicle, a supersonic combustion ramjet (scramjet) en- 
gine was proposed to provide the power. Inside this 
scramjet engine, compressible mixing layers are impor- 
tant phenomena. The performance of the engine will 
depend on the supersonic mixing and the flame holding 
of shear layers. 

The behavior of incompressible mixing layers has 
been studied extensively. However, additional study 
is required to understand the compressibility effects of 
free shear layers at high speeds. Figure 1 shows a sketch 
of a typical free shear layer. Two streams at different 
temperatures, densities, and Mach numbers merge to- 
gether to form a free shear layer. Various combinations 
of flow conditions of high-speed streams and low-speed 
streams allow for the systematic study 
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of compressible shear layers. In this paper, we report 
the incorporation of a k — € model with compressibility 
effects to a two-dimensional flow equations solver for 
the simulation of compressible shear layers. First, we 
point out the derivation procedure of the compressible 
k and c equations. Unlike the procedure in the incom- 
pressible equations, both Favre and Reynolds averaging 
procedures 1 are performed in deriving the flow and tur- 
bulence equations. Particularly, the additional terms in 
the Navier-Stokes equations due to the averaging pro- 
cedure are illustrated. These terms are often orruued 
in CFD practices. The k and e equations are presented 
in vector form for the convenience of illustrating the 
numerical method. 

The lower-upper (LU) scheme developed by Yoon 
and Jameson 2 is adopted in this work. This method has 
proven very efficient for large systems of equations. 3 For 
completeness, a brief account of the numerical method 
is presented in this paper. The newly developed solver 
then is applied to simulate compressible free shear lay- 
ers with five different convective Mach numbers from 
0.05 to 1.48. One of five test conditions is a replica 
of that reported by Goebel and Dutton. 4 The most 
important feature of the compressible free shear lay- 
ers one want to demonstrate in the calculations is the 
decrease of the shear layer thickness with increasing 
Mach number. 5 Favorable comparison were found be- 
tween the experimental data and the calculated results. 


THEORETICAL MODEL 

In deriving the compressible flow equations of fully 
turbulent flows, all the flow properties are Favre aver- 
aged (mass weighted averaged) except the density, p, 
and pressure, p. The definition of the Favre average is 

<t> = <)>pI~p - 
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(1) 


where <f> is any flow property Thus flow variables are 
decomposed in the following fashion, 


<t> — 4> + <i> n ■ 


( 2 ) 


On the other hand, conventional Reynolds’ average are 
used for the pressure and the density. According to 
the definition of Favre averaging, the following relations 
exist: 


r = - 


pW 

p 


P<F = 0, 


(3) 


4>* = o. 


In doing so, all the terms associated with the density 
fluctuation, e.g., p f u * , in the Reynolds’ averaged equa- 
tions were eliminated in the Favre ’s averaged equations. 
The resulting flow equations are much simpler com- 
pared to the equations derived by the Reynolds aver- 
aging procedures. 

Written in a strong conservative form, the turbu- 
lent, compressible, flow equations can be expressed as 
follows: 


^ + |( E -E„)+| ; ( F -F,) = H. (4) 

Here x and y are Cartesian coordinates, Q is the depen- 
dent variable, E and F are the convective flux vectors 
and E v and F v are the viscous flux vectors. The equa- 
tions are similar to the laminar equations. However, 
all variables in the equations are the averaged vari- 
ables, and the transport properties are the effective, 
i.e., laminar plus turbulent, properties. Here, we want 
to point out that the viscosity multiplied by the dilata- 
tions! terms in the normal stresses is still the laminar 
viscosity as illustrated in the following relations: 


T X x 

T yy 


.du 2 f du dv\ 
dv 2 f du dv\ 


(5) 


The vector H on the right hand side of the flow 
equations, Eq. 4, represents the additional terms intro- 
duced by the averaging procedure. The vector H can 
be expressed as 


H = 


/ 0 \ 

3 dx 

ai [(/* + £)fl] + ^ [O' + £)fj] 

+2«£(?*) + 25£(p*) + G / 


( 6 ) 


where k is the turbulent kinetic energy and is defined 
as 

k = ^p(u"u +v"v"), ( 7 ) 

2p 

G is the generation of the turbulent kinetic energy and 
is defined in the following section. 

The equations of turbulent kinetic energy, k , and 
energy dissipation rate, c t are derived by the manipula- 
tion of the flow equations and the averaging procedure. 
The derived equations of k and c can not be solved di- 
rectly due to the closure problem. Modelling of certain 
terms in the equations is necessary to make the govern- 
ing equations well posed. The details of the modelling 
is beyond the scope of this work. Here, only the final 
form of the equations are presented. The k and c equa- 
tions in the Cartesian coordinate system can be cast 
into the vector form: . 


dQfcc dEjfce dF ke 

dt + dx dy 


^E»/Jbe , dF vic t c 

~dT + d y + ’ 


( 8 ) 


where 



V (/*+£)$!’ 

( G — + atMf)\ 

V (CiG-CaAOi ) ’ 


E v t e — 
Fvke ~ 


(9) 


where G is the generation term of the turbulent kinetic 
energy and can be expressed as: 


G ~ p t 



2 du k 

3 dxt . 


2 

3 pt ar t 


( 10 ) 


where the subscripts :, j y and k follows the convention 
of Cartesian tensor. 

The eddy viscosity p t is derived in terms of appro- 
priate length and velocity scales. For the Jb — e turbu- 
lence model, the length scale and the velocity scale of 
turbulent fluctuations are taken as 



it 1 ' 2 . 


(ii) 
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These relations allow the eddy viscosity, p tt to be mod- 
elled as: 

_k 2 

Vt = C^p — . ( 12 ) 

The constants used in the k - c model are the standard 
Jones and Launder *s values: 6 C \ - 0.09, C x = 1.44, 
O 2 — 1.92, o'* = 1.0, and <x t = 1.3. These constants 
were never altered during the course of this work. 

In order to accommodate the compressibility ef- 
fect, the dissipation term pc in the source term of the k 
equation is multiplied by a correction factor (1 -f orAf 2 ). 
Here M t is the local turbulence Mach number defined 
as M t = y/k/a where a is the local speed of sound. The 
constant a in the term is taken as unity. This model 
is developed by Sarkar et al. 7 The physical meaning of 
the term is that for high turbulence Mach number ( M t ) 
flows, the dissipation of the turbulence kinetic energy 
is enhanced by a factor of crA/ 2 . For free shear layers at 
high convective Mach number, the turbulence intensity 
is greatly reduced due to compressibility. 

In calculating the turbulent free shear layers, the 
inlet boundary conditions for mean velocities' and tem- 
perature are specified based on the hyperbolic tangent 
profile with specified initial shear layer thickness. The 
hyperbolic tangent profile is an approximation of the 
self-similar solution for fully developed turbulent free 
shear layers. The inlet transverse velocities are set to 
be zero. The turbulent kinetic energy and dissipation 
are specified according to the local equilibrium assump- 
tion and a algebraic turbulence model: 8 


- -S/2 d “ 


plm dy 


(13) 


where l m = 0 .1256 and the shear layer thickness, 6, is 
based on the distance between the two transverse lo- 
cations where u = tii - O.lAu and u = u 2 + O.lAtt. 
The dissipation can be related to the local length scale 
which is specified based on the local shear layer thick- 
ness: 

- ** 

€ = c <t ( 14 ) 

where C t = 1.23. Using Eqs. (13) and (14), and the 
equation for the eddy viscosity, i.e, Eq. (12), k and 
€ can be readily obtained for the upstream boundary 
conditions. 


Numerical Method 

The numerical solution of Eqs. (4) and (8) is 
performed in a general, body-fitted coordinate system, 
(£»*?)• For the purpose of discussion, we will concen- 
trate on Eq. (4). However, the procedure is equally 


applicable to Eq. (8). Coordinate transformation of 
Eq. (4) results in: 


5Q d / - ~\ d / ~ - \ 

57 + Ti ( E - E 0 + Tri ( f - F 0 = H ( 15 > 

where 

Q = hQ 

E = h(f x E + ^F) 

F = h( Vx E + r? y F) 

E„ = -K y F„) (16) 

F„ = h(ij z E„ -f r? y F v ) 

H = hH 

in which h is the cell volume. 

The transformed equation, Eq. (15), is solved us- 
ing a time-marching, LU scheme. The LU scheme can 
be obtained by approximately factorizing the left-hand- 
side (LHS) of the equation. In time-marching form, the 
implicit upwind difference scheme of Eq. (15) can be 
written as 


(I + At(D+A~ + D+B~ - D+ 

A+ + D~ B + )]AQ = At RHS' ^ 

In Eq. (17), At is the time-step. Backward-difference 
operators are denoted by D~ and Z)“, and forward- 
difference operators by D+ and D+. The flu* Jaco- 

bians, A + , B + , A~, and B~ are constructed such that 
the eigenvalues of matrices are non negative and 
those of matrices are non positive. The matrix D 
is the Jacobian matrix of the source term. 

The LHS matrix of Eq. (21) is usually too large for 
direct inversion. An approximate-factorization proce- 
dure is implemented which results in the following LU 
scheme : 


I+A< / 

f 

D 7A+ 

+ D~ B+ - 

A~ 

B“ - V 

D 

L \ 
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\ 

7 

A£ 

A, )\ 

= A<R 






AQ 


where the grid spacing in the general coordinate, A£ 
and At) are usually taken to be one. R represents the 
residual of each LU time marching step. In calculating 
the residual, R, both the inviscid and viscous terms are 
discretized using the central-difference approach: 

R = £>t( E v - E) + D V (F V -F) + H (19) 
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where D $ and D n are the central difference operators. 

Equation (18) is the generic form for the LU 
scheme. Its derivation can be found in Ref. 3 and 
will not be repeated here. This LU scheme requires 
inversion of the matrix, 


is applied in this work to enhance the numerical stabil- 
ity. In the k equation, c has been replaced by k(c/k) 
where c/k is treated as a constant. A similar method is 
also applied to the source term of the c equation. The 
Jacobian matrix obtained is: 


I + A< (A + - A - + B+ - B- - 6)] (20) 


/_x (1 + qM 2) o \ 

V o -C,iJ 


(27) 


for the L operator and 

[i + A* (A+ - A" + B + -B“ — &/£>)] (21) 

for the U operator. 

Up to this point, no definition has been made to 
the exact form of the split flux Jacobians. Yoon and 
Jameson 2 proposed that the split flux Jacobians are 
defined as 

A+ = 0.5(A + 7A I) 

A" = 0.5(A- 7 :I) 

- ( 22 ) 

B+ = 0.5(B + 7gl) 

B- = 0.5(B - 7*1) 

where 7^ and 7^ are greater than the spectral radii of 
the associated flux Jacobians : 


7a > ™x(|A A |) 

7b > m ax( | Aq 1) 

The purpose of constructing split flux- Jacobians by Eq. 
(22) is to make the matrices in Eqs. (20) and (21) diag- 
onal for efficient inversion. Apparently, the eigenvalues 
of the split flux- Jacobians are not the characteristics 
speeds of the flow. 

In solving the k and c equations, the aforemen- 
tioned numerical method, i.e., the LU scheme on the 
left hand side and central differencing on the right hand 
side, is used. The solution procedure of the whole equa- 
tion set is decoupled into flow solver and turbulence 
solver. Thus, the turbulence solver stands alone and 
can be easily turned on or off. This arrangement does 
not affect the overall numerical stability due to the fact 
that the feedback from k and c equations to the flow 
equations depends on the turbulent transport proper- 
ties only. Thus, it is more efficient and convenient to 
separate the solution procedure into two parts. 

The source terms of the k and c equations demand 
special treatment. In linearizing the source terms for 
the numerical method, the Jacobian matrix is obtained 
through the derivative of the source terms with respect 
to the dependent variables, i.e., pk and pc. Following 
the usual practice, the form of the source terms guar- 
antee a 2 x 2 full matrix for the Jacobian matrix. How- 
ever, special treatment in deriving the Jacobian matrix 


Note that off-diagonal terms are eliminated and the di- 
agonal terms are always negative. Thus, the implicit 
part of the source terms of k and c equations behaves 
like a sink which always stabilize the numerical scheme. 


Results and Discussions 

Bogdanoff 5 introduced the convective Mach num- 
ber as a parameter that collapses the growth rate date 
of plane shear layers. The convective Mach number is 
defined as 


Ui - Uc = Uc - U2 
ci c 2 


(28) 


where U\ and t/ 2 are the freestream velocities, and 
and c 2 are the freestream sound speeds. U c is the con- 
vective velocity of large structures and is defined as 


U c = 


U ic 2 -f U2C1 

Ci + c 2 


(29) 


Table 1 shows the five test conditions of the simulated 
compressible free shear layers. The range of the convec- 
tive Mach numbers, M c , is from 0.05 to 1.48. The last 
two rows show the spreading rate and the ratio of the 
compressible spreading rate to incompressible spread- 
ing rate in which 6 is the shear layer thickness and the 
superscript ' represents the derivative of 6 with respect 
to the streamwise distance. As indicated in the last 
row of the table, the ratios of the spreading rate of 
the free shear layers decreases as the convective Mach 
number increases. The test conditions of Case 3 in the 
table are the same as in the experiments reported by 
Goebel and Dutton. 4 In the rest of the section, we will 
first show the direct comparison between the simulated 
results and the experimental data for Case 3. Then, 
detailed numerical solutions of Case 3 in terms of ve- 
locity, turbulent kinetic energy, turbulence dissipation 
rate, Reynolds stress, and eddy viscosity are presented 
in a coalesced fashion. Finally, the solutions of five 
cases are compared to each other in the figures of ki- 
netic energy, Reynolds stress, and ratio of spreading 
rates. 

The solutions of Case 3 are examined in detail. 
Figure 2 shows the development of the free shear layer. 
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Note that x and y axes are not on a 1:1 ratio for the 
convenience of illustration. The definition of the shear 
layer boundaries is the same that of the shear layer 
thickness. In Fig. 2, the boundaries of the shear layer 
corresponding to 10 - 90% are drawn. Circles are the 
experimental data of Dutton et al. The calculated re- 
sults agree well with the experimental data. After the 
developing region, the boundary of the shear layer is 
almost linear. Incidentally, Figs. 3a and 3b show the 
numerical convergence trends of the flow and it — c 
equations. In about 2500 iterations, the residuals drop 
about 12 orders of magnitude and reach the machine 
accuracy. 

Figure 4 shows the Mach number profile at various 
axial locations. The velocity gradient in the transverse 
direction decreases as the flow goes downstream. Figure 
5 is the coalesced version of Fig. 4. The nondimension- 
alized y coordinate defined as (t/— y c )/<5 is used, where 
y c is the transverse location at the center of the shear 
layer, and S is the shear layer thickness. Note that the 
upstream boundary condition of velocities is prescribed 
according to a hyperbolic tangent curve which is an ap- 
proximation of self-similar solution of free shear layers. 
According to Fig. 5, this self similarity of velocity pro- 
files never fail as the flow goes downstream. 

Figure 6 shows coalesced turbulence kinetic ener- 
gies at various locations. The turbulence kinetic en- 
ergy is nondimensionalized by AU 2 . This figure clearly 
shows that after the first three stations, i.e., about 150 
mm, the turbulence kinetic energy retains self similar- 
ity. Thus, the developing region for the turbulence is 
about 150 mm. Figure 7 shows the turbulence dissipa- 
tion profiles. As flow goes downstream, the peak values 
of turbulence dissipation at each stations decrease, and 
the turbulence dissipation never reach a fully developed 
condition. However, if the € value is nondimensional- 
ized by A U 3 /6 the coalesced profiles appear as shown 
in Fig. 7. A similar behavior is observed for the eddy 
viscosity profiles (Fig. 8). Figure 9 shows the nondi- 
mensionalized Reynolds stress profiles. Again, the tur- 
bulence Reynolds stress becomes the fully developed at 
about 150 mm downstream of the splitter plate. 

Figure 10 shows the comparison between the pre- 
dicted fully developed Reynolds stress and the experi- 
mental data reported by Dutton et al. 5 The predicted 
solution underestimated the peak value of the Reynolds 
stress profile by 6 - 8%; however, the overall trend 
of the predicted results is correct. Many factors con- 
tribute to the discrepancy between the predicted result 
and experimental data. Among them, the upstream 
boundary conditions are simplified in the solution pro- 
cedure, i.e., no effort was made to simulate two bound- 
ary layers merging at the tip of the splitter plate. This 
could offset the solution in the developing region and 


shift the fully developed solutions. 

Figure 11 shows the comparison of the turbulent 
kinetic energy between the five cases. Again, the turbu- 
lent kinetic energy is normalized by the square of the ve- 
locity difference of the two streams. It is clear that the 
normalized turbulence intensity decreases with increas- 
ing convective Mach number. A similar situation is 
observed in Fig. 12, the normalized Reynolds stress de- 
creases with increasing convective Mach number. Fig- 
ure 13 shows the distribution of the ratios of spread- 
ing rates for the five cases compared to experimental 
data. The ratio of the spreading rates decreases from 
about unity to 0.45 as the convective mach number in- 
creases from about 0. to 1.45. Both the experimental 
data and the simulated results show the spreading rate 
ratio reaches an asymptotic value after the convective 
mach number exceeds unity. Simular phenomenon can 
be seen in Figs. 11 and 12. Both figures show that 
the normalized turbulence kinetic energy and Reynolds 
stress reach asymptotic values as the convective Mach 
number increases. 


CONCLUDING REMARKS 

In this paper, we report the incorporation of a com- 
pressible k — € model with a two-dimensional Navier 
Stokes solver to study compressible free shear layers. 
Sarkar’s modelling is adopted to simulate the compress- 
ibility effects of the k and e equations. The model en- 
hances the turbulence dissipation rate of flows at high 
speeds. In deriving the governing equations, the Favre 
averaging procedure for the fully trubulent flow equa- 
tions is elaborated. The equation sets are presented in 
the vector form for the convenience of the discussion of 
the numerical method. Yoon and Jameson’s LU scheme 
is used to solve the equation sets. Details of boundary 
conditions and the treatment of the source terms of the 
k and e equations are discussed. Then the program is 
applied to simulate compressible free shear layers with 
five different convective Mach numbers. The decrease 
of the spreading rate of the shear layers with increasing 
Mach number is observed in the calculated results. Re- 
sults also show favorable comparison with Goebel and 
Dutton’s experimental data. 
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Fig. 1: Schematic of free shear layer. 


Table 1: Test conditions of the calculated shear layer. 


CASE 

1 

2 

3 

4 

5 

m,,m 2 

1.2, 1.1 

2.0, 1.4 

1.91, 1.37 

2.5, 1.1 

6.1, 3.15 

M c 

0.05 

0.31 

0.45 

0.70 

1 .48 

u lt u 2 

[m] 

419, 384 

676, 465 

702, 404 

865, 380 

3461, 1786 


300, 300 

275, 275 

334, 215 

295, 295 

800, 800 

(K?g/m^3] 

0.64,0.64 

0. 7,0.7 

0.57,0.89 

0.64,0.64 

0.24,0.24 

P [atm] 

0.55 

0.55 

0.55 

0.55 

0.55 

8’ 

0.007 

0.021 

0.027 

0.035 

0.024 

878’ j 

1.01 

0.7 

0.54 

0.47 

0.40 
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Log (Iteration error) Lofldteritlon error) i Vertical Location 


15 mm 




Fsg. 4 : The velocity profiles of the free shear layer at dif- 
ferent axial locations. 


Fig. 2: The boundaries of the free shear layer (10 90%). 
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Fig. if: Coalesced velocity 
tions. 
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profiles of different axial loca- 


Fig. 3: Convergence trends of the flow and turbulence 
equations. 



















F'g- 11: Comparison of turbulent kinetic energy with five 
different convective Mach numbers. 
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Papamoschou and Roshko [9] 
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Langley’s data [7] 
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Fig. 13: Comparison of ratios of spreading rate between ex- 
perimental data and predicted results. 


Fig. 12: Comparison of Reynolds stresses with five different 
convective Mach numbers. 
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